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Abstract. We study the distinction and quantification of chaotic and regular motion 
in a time-dependent Hamiltonian barred galaxy model. Recently, a strong correlation 
was found between the strength of the bar and the presence of chaotic motion in this 
system, as models with relatively strong bars were shown to exhibit stronger chaotic 
behavior compared to those having a weaker bar component. Here, we attempt to 
further explore this connection by studying the interplay between chaotic and regular 
behavior of star orbits when the parameters of the model evolve in time. This happens 
for example when one introduces linear time dependence in the mass parameters of 
the model to mimic, in some general sense, the effect of self-consistent interactions of 
the actual N-body problem. We thus observe, in this simple time-dependent model 
also, that the increase of the bar's mass leads to an increase of the system's chaoticity. 
We propose a new way of using the Generalized Alignment Index (GALI) method 
as a reliable criterion to estimate the relative fraction of chaotic vs. regular orbits 
in such time-dependent potentials, which proves to be much more efficient than the 
computation of Lyapunov exponents. In particular, GALI is able to capture subtle 
changes in the nature of an orbit (or ensemble of orbits) even for relatively small time 
intervals, which makes it ideal for detecting dynamical transitions in time-dependent 
systems. 
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1. Introduction 

The study of chaotic and regular properties of the motion in Hamiltonian systems 
constitutes a vast area of research in the field of nonlinear dynamics. Since the early 
1960's, several methods and tools for the fast and accurate detection of the nature of 
orbits have been proposed and applied to this end in a great number of publications. 
One may refer e.g. to the pioneering paper by Henon and Heiles pQ, where the Poincare 
Surface of Section (PSS) (2j section 1.2b] was used to reveal the chaotic properties of 
a non-integrable 2 degree of freedom (2 dof ) Hamiltonian system. Of great importance 
in this direction was also the algorithm proposed by Benettin and co-workers [31 HI [5] 
regarding the computation of the full spectrum of Lyapunov exponents (LEs) associated 
with the time evolution of a single deviation vector from a reference orbit, which applies 
to dynamical systems of arbitrary dimension. More recently, other related methods 
have been proposed in the literature, like the "Fast Lyapunov Indicator" [61 [7] and 
the "Mean Exponential Growth of Nearby Orbits" (MEGNO) 0,19], while there have 
also been approaches focusing on the time series constructed by the coordinates of each 
orbit, like the "Frequency Map Analysis" [HS CEIl [32] and the "0-1" test [331 El EH]. 
Interesting accounts of these methods can be found in [16], as well as in a more recent 
review paper [5]. 

A novel, very efficient method based on the evolution of k > 2 initially linearly 
independent deviation vectors is provided by the so-called "Generalized ALignment 
Indices" (GALI or GALI& spectrum) introduced in [T7] as a generalization of the 
"Smaller ALingment Index" (SALI) [HI [191 [20]. The major advantage of the GALI 
method is that it follows the evolution of two or more deviation vectors and is thus 
able to extract more information about the complexity of the motion, yielding i.e. the 
dimensionality of the invariant torus on which a regular orbit lies and predicting faster 
the chaotic nature of trajectories [211 E21 [23]. To date, the GALI and the SALI indices 
have been successfully applied to a wide variety of autonomous (i.e. explicitly time- 
independent) conservative flows and maps (see e.g. [24| [25] [26] [27[ [28] [29] [30] [8T[ [32] 
[MlElESllMlETlEHlEgiOlIlT]). A concise review of the theory and applications of 
both the SALI and GALI methods can be found in [4*2l Chapter 5]. 

The motivation of the current work is twofold: First, we wish to investigate the 
different dynamical properties of a non-autonomous galactic potential, whose time- 
dependence could mimic certain realistic general trends arising in barred N-body galaxy 
simulations. Our second main goal is to explore the advantages of the GALI method, 
over the more traditional LEs, in detecting dynamical transitions in Hamiltonian 
systems, whose equations of motion are explicitly time-dependent. 

There are, of course, several studies of time-dependent (TD) galactic and 
cosmological models in the literature, which use different tools to identify the chaotic 
vs. regular nature of orbits. Defining the orbital complexity n(k), of an orbital segment 
as the number of frequencies in its discrete Fourier spectrum that contain a fc-fraction 
of its total power [43] . one may compare n(k) with the short-time evolution of the LEs 
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for TD models [44]. In [15] the case of a cosmological model is discussed, where orbits 
may experience regular and/or chaotic motion during their time evolution, while in [16] 
the effects of a black hole, friction, noise and periodic driving are studied on a triaxial 
elliptic galaxy model, in which a type of transient chaos was found caused by a damped, 
oscillatory component [U_l EE] . Finally, in [19] the so-called "pattern method" was used 
to study a Henon-Heiles potential to which an exponential function of time is added, 
while the dynamics of some simple TD galactic models was investigated in [501 EI]- 

We recall here that in conservative systems the asymptotic nature of an orbit may 
be either periodic, quasiperiodic or chaotic. In the latter case, however, it may take a 
very long time before one can safely claim that a "final" state is reached, depending on 
the local dynamical properties, which may be characterized by "strong or weak chaos" . 
Here we explore the advantages of the GALI method and compare its predictions with 
what one finds using more traditional methods like the computation of the maximum 
Lyapunov exponent (MLE). 

In particular, we focus our attention on the dynamics of a barred galaxy model 
containing a disc and a bulge component, which is a widely accepted model for real 
barred galaxies. In the spirit of a mean field approach, we consider the motion of 
stars (represented by point particles) in this potential. The richness of the dynamics 
of the time-independent (TI) version of this model has been extensively studied in 
terms of: (a) the detection of periodic orbits and the analysis of their stability (see 
e.g. [531 EU ESI ESI E3 EH EH]), (b) the estimation of the relative fraction of chaotic 
vs. regular orbits [321 EOj [37] , and (c) the statistical distributions of orbital coordinates 
described by g-Gaussian distribution functions [61]. 

Here, we extend the analysis by considering a TD version of this model. More 
specifically, we allow some mass parameters of the potential to vary linearly as functions 
of time. As expected, whether we study a 2 dof or 3 dof version of the model, these 
variations can change the stability properties of periodic orbits, "dissolve" islands 
of regular motion and alter the structure of phase space in very complicated ways. 
Recently, it was found in the TI case that the relative fraction of chaotic orbits grows 
as the bar's strength increases [37]. A question therefore arises, whether a similar 
correlation holds in the TD model, in the presence of "realistic" trends, which permit 
the mean field potential to vary in a way that is compatible with self-consistent N-body 
simulations regarding several components of the system. 

Clearly, the analysis of the full N-body problem describes much better the galactic 
evolution and captures in great detail the different stellar structures present in the 
dynamics. However, there are serious difficulties and limitations when one tries to 
apply dynamical chaos detectors to such "realistic" many-particle systems due to the 
lack of sufficient orbital information during the time evolution. For this reason, many 
researchers prefer to use mean field potentials that are "frozen" in time and study the 
properties at specific snapshots of the simulations [621 [63] . 

Keeping in mind that a barred galaxy experiences several dynamical transitions 
in different epochs that cannot be easily incorporated in our TD mean field potential, 
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we shall proceed by making some helpful assumptions in an attempt to understand the 
behavior of such widely used chaos detectors as the GALI and the MLE. Thus, we will 
treat here two very general dynamical trends known to occur in barred galaxies: In the 
first scenario the mass of the bar component grows linearly in time (at the expense of 
the disc mass). This increase may be caused by an exchange of angular momentum 
with the disc (outer parts gain momentum from the inner parts), as has already been 
observed in N-body simulations (see e.g. [64, 65J). The fundamental trend in this case 
is that bars generally grow stronger in time, become more elongated and massive and 
eventually slow down. We will also consider the inverse scenario, where the bar gets 
weaker making the disc more massive as time evolves (see e.g. [66], IBTj). 

The paper is organized as follows: In section [2] we present the TD barred galaxy 
model used in our study, while section [3] is devoted to the description of the numerical 
methods employed for the computation of the MLE and the GALIs. Section H] contains 
the main numerical results of the paper. A detailed investigation of the dynamics of 
particular orbits in a 2 dof version of our model is performed in section 14.11 while orbits 
of the full, 3 dof model, are considered in section 14.21 A global investigation of the 
dynamics of our TD galactic model is given in section fl~3l Finally, in section [5] the main 
conclusions of our work are presented. 



2. The model potential 

Let us consider the following TD 3 dof Hamiltonian function which determines the 
motion of a star in a 3 dimensional rotating barred galaxy: 

H = o(P* + Py + Pi) + V ( x > V' z > *) -tt b {xp y -yp x ). (1) 

The bar rotates around its z-axis (short axis), while the x direction is along the major 
axis and the y along the intermediate axis of the bar. The p x , p y and p z are the 
canonically conjugate momenta, V is the potential, fib represents the pattern speed of 
the bar and H is the total energy of the orbit in the rotating frame of reference (equal 
to the Jacobi constant in the TI case). 

The corresponding equations of motion are: 

x = p x + Q b y, 

y = Py~ flbX, 
Z = Pz, 

dV 

P* = ~-^ + nb Pv ( 2 ) 

dv 

Py = MbPx, 

oy 
dV 

P* = o ' 

OZ 

while the equations governing the evolution of a deviation vector w = 
(5x, 5y, 5z, 5p x , 5p y , Sp z ) needed for the calculation of the MLE and the GALIs, are 
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given by the variational equations: 
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6 

dzdy 



d 2 V 

dxTz 5z + ^ ( 3 ) 
d 2 V 

— — -6z - VL b 5p x , 
dydz 

d 2 v 5z 

dzdz 

The potential V of our model consists of three components: 

A triaxial Ferrers bar [52], the density p(x,y, z) of which is given by: 

i \ I Pc(l-^ 2 ) 2 if m < 1, ,.s 
p(x,y,z) = < p v ; „ (4) 

HK ' y ' ; 1 if m> 1, K J 

where p c = |f^ gj ^ is the central density, Ms{t) is the mass of the bar which 
changes in time, and m 2 = fj + ^j + fj, a > b > c > 0, with a, b and c being the 
semi-axes of the ellipsoidal bar. The corresponding potential is: 

V B = -nGabc^ J™ ^(1 - m 2 {u))\ (5) 
where G is the gravitational constant (set equal to unity here), m 2 {u) = -^r^ + 

2 2 

j^— + -^t^i A 2 (u) = (a 2 + u)(b 2 + u)(c 2 + u), and A is the unique positive solution 
of m 2 (A) = 1, outside of the bar (m > 1), while A = inside the bar. The analytical 



expression of the corresponding forces are given in 

(b) A bulge, modeled by a Plummer sphere [68] whose potential is: 

Vs = GM * ( 6 ) 

^ x 2 + y 2 + z 2 + el 

where e s is the scale-length of the bulge and Ms is its (constant) mass. 

(c) A disc, represented by the Miyamoto-Nagai potential 



V D = - , GM ° (t) (7) 



'x 2 + y 2 + (A + y/z 2 + B 2 ) 2 

where A and B are its horizontal and vertical scale-lengths, and the mass of the 
disc Mo(t) changes in time so that the total mass of the system is kept constant. 

The model's parameters have the following constant values: G — 1, f2b=0. 054 (54 
km ■ sec -1 ■ kpc~ x \ a— 6, b =1.5, c=0.6, A=3, B=l, Ms=0.08, while the initial values of 
the bar and disc masses are Mb(0)=0.1 and Md(0)=0.82, respectively. The units used 
are: 1 kpc (length), 1000 km- sec -1 (velocity), 1 Myr (time), 2 x 10 11 Mq (mass). The 
total mass Ms + Mn(t) + Ms(t) is set equal to 1 and since the bulge's mass Ms is kept 
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constant, the disc's mass M£>{t) is varied as Md(£) = 1 — (Mg + Mg(t)). The rate of 
the mass variation of the bar is chosen to be linear according to the law: 



where the proportionality constant is a > or a < respectively, if the mass of the bar 
increases or decreases in time. 

In order to measure the time variation of the bar's strength we calculate the quantity 



which estimates the relative strength of the non-axisymmetric forces. In the above 
expression, $ is the potential on the symmetry plane z = expressed in polar 
coordinates (r,8), <3> is its axisymmetric part, while the maximum in the first term 
on the right hand side of (Q is calculated over all values of the azimuthal angle 9. The 
maximum value of Qt(r) over all radii shorter than the bar extent, termed Qf>, can be 
used as a measure of the bar's strength. 

It is clear, of course, that the variation of the bar's strength modifies the values of 
several parameters and yields richer information about the dynamics of a self-consistent 
model. N-body simulations show that in general, variations of the bar's mass also 
change the mass ratios of the model's components, the bar's shape and the pattern 
speed of the galaxy. Hence, if one wishes to use a mean field potential to "mimic" a self- 
consistent model as accurately as possible, one should allow for all the parameters that 
describe the bar (together with all other axisymmetric components) to depend on time, 
assuming that the laws of such dependence were explicitly known. In our case, however, 
we adopt a simpler approach and vary only the masses of the bar and the disc, as a 
first step towards investigating such models when time dependent parameters are taken 
into account. Thus, we do not pretend to be able to reproduce the exact dynamical 
evolution of a realistic galactic simulation. Rather, we wish to understand the effects 
of time dependence on the general features of barred galaxy models and compare the 
efficiency of indicators like the GALIs and the MLE in helping us unravel the secrets of 
the dynamics in such problems. 

3. Computational methods 

In order to estimate the value of the MLE, Ai, of a particular orbit we follow the 
evolution of the orbit and a deviation vector w from it, by numerically solving the set 
of equations (j2J) and (131) respectively. For this task we use a Runge-Kutta method of 
order 4 with a sufficiently small time step (typically of the order of r w 10~ 2 ), which 
guarantees the accuracy of our computations (i.e. the use of the half time step does not 
practically change our results). 

The ordinary differential equations (ODEs) (J2]) can be solved independently from 
equations (EJ). On the other hand, the latter set of equations, governing the evolution 



M B {t) = M fl (t = 0) + at, 



(8) 



[701 HQ: 




(9) 
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of a deviation vector, has to be solved simultaneously with equations because the 
second derivatives of the potential V, appearing in the right hand side of ([3]), depend 
explicitly on the solution of (T5]). Note that (T2]) constitutes a non-autonomous set of 
ODEs because the derivatives of V depend explicitly on time. Although one could 
transform (T5]) (and consequently ([3])) to an equivalent autonomous system of ODEs by 
considering time t as an additional coordinate (see e.g. [21 section 1.2b]), this approach 
is not particularly helpful, and is better to be avoided [72] . 

So, in order to compute the MLE and the GALIs we numerically solve the time- 
dependent set of ODEs ©and ([3D . Then, according to [731 H2 El H the MLE A x is 
defined as: 



Ai = lim oi(i), (10) 



t— >oo 



where: 



is the so-called "finite time MLE", with ||w(0)|| and ||w(£)|| being the Euclidean norm 
of the deviation vector at times t = and t > respectively. A detailed description of 
the numerical algorithm used for the evaluation of the MLE can be found in [5]. 

This computation can be used to distinguish between regular and chaotic orbits, 
since Ci(t) tends to zero (following a power law oc t^ 1 ) in the former case, and converges 
to a positive value in the latter. But Hamiltonian (TjQ) is TD, which means that its 
orbits could change their dynamical behavior from regular to chaotic and vice versa, 
over different time intervals of their evolution. In such cases, the computation of the 
MLE (fit)]) might not be able to identify the various dynamical phases of the orbits, since 
by definition it characterizes the asymptotic behavior of an orbit. 

In order to avoid such problems in our study, we also turn to the use of the GALI 
method of chaos detection [T7J. The GALI index of order k (GALI&) is determined 
through the evolution of 2 < k < N initially linearly independent deviation vectors 
w.j(0), i = 1,2, ...,k, with N denoting the dimensionality of the phase space of our 
system. Thus, apart from solving (J2j) , which determines the evolution of an orbit, we 
have to simultaneously solve ([3]) for each one of the k deviation vectors. Then, according 
to [17] . GALIfc is defined as the volume of the fc-parallelogram having as edges the k 
unit deviation vectors Wj(t) = Wi(t)/||wj(t)||, i = 1,2, ...,k. It can be shown, that this 
volume is equal to the norm of the wedge product (denoted by A) of these vectors: 

GALI fe (t) =|| wi(t) A w 2 (t) A ... A w p (t) || . (12) 

We note that in the above equation the k deviation vectors are normalized but their 
directions are kept intact. In practice, we apply a numerical method for calculating 
this norm, which is based on the singular value decomposition of an appropriate matrix 
[751 EQ. 

The behavior of GALI& for regular and chaotic orbits was theoretically studied in 
[17] [21] . where it was shown that all GALI&(t) tend exponentially to zero for chaotic 
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orbits, with an exponent that depends on the first k LEs of the orbit. In the case of 
regular orbits, GALI^ remains practically constant and positive if k is smaller or equal 
to the dimensionality of the torus on which the motion occurs, otherwise, it decreases 
to zero following a power law decay. 

In order to use GALI^ to capture the dynamical changes of orbits in TD systems we 
apply the following procedure: Whenever GALI^ reaches very small values (i.e. GALIj, < 
10 -8 ) we reinitialize its computation by taking again k new random orthonormal 
deviation vectors, which means that we set again GALL; = 1. Then we let these vectors 
evolve under the current dynamics. An exponential decrease of GALL indicates chaotic 
behavior. Thus, the time tj needed for GALL to become less than 10~ 8 can be used to 
identify epochs where the orbit is chaotic or regular. 

4. Numerical results 

4-1. A typical orbit of the 2 dof model 

Let us first start with the simple example of an orbit whose motion is restricted in the 
2 dimensional (2D) (x,y) space (or 4D phase space), of a 2 dof reduced version of the 
full 3 dof model, where z,p z are set equal to zero at t = 0, and remain zero at all 
times. Furthermore, we shall assume that the bar component gets stronger in time, 
according to (JSj) , in agreement with the general trend in the dynamical evolution of 
barred galaxies (e.g. [6J1 [65]). In general, the fraction of the chaotic component here 
is expected to increase as the bar gains mass and becomes stronger [37J. However, this 
does not necessarily imply that the nature of all orbits remains unchanged in time; 
regular orbits can become chaotic as time goes on, and vice versa. For this system one 
can construct, at specific times, 2D phase portraits (similar to the PSS of the TI case) 
to visualize the dynamics and follow the changes that orbits undergo. 

Assuming that all other components except M#(i) and Mo(t) remain constant, we 
vary Mb from Ms(t = 0) = 0.1 to MB{tjinai = 20000) = 0.2, yielding a proportionality 
constant a = 5 x 10~ 6 in ([8]). The initial energy value for these parameters is set to be 
HAito) = E^ito) = —0.2570. The strength of the bar, measured by the Qb parameter is 
initially Qb(t ) = 0.425 and becomes Qb(tfi na [) = 0.6732. 

As a test case, we pick a specific orbit, which we call orbit A, with initial condition 
(x,y,p x ,p y ) = (0.0,1.0,0.16531,0.0). This orbit is initially located inside an island of 
stability in the system's phase space and is thus expected to remain regular, at least for 
some time to come. 

In the left column figure [1] we show the projections of orbit A on the (x, y) plane, 
while on the right column we plot its intersection points with a PSS defined by x = 0, 
Px > (black points) for five successive time intervals, each having a duration of 2500 
time units. In every panel of the right column we plot in gray the PSS which corresponds 
to the central time of each time interval, i.e. t = 1250, t = 3750, t = 6250, t = 8750 and 
t = 11250, respectively from top to bottom. Orbit A is initially located inside the right 
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Figure 1. Left column: projections of orbit A on the (x, y) configuration space for 
successive time intervals with 2500 time units length: < t < 2500 (first row, interval 
1), 2500 < t < 5000 (second row, interval II), 5000 < t < 7500 (third row, interval III), 
7500 < t < 10000 (fourth row, interval IV) and 10000 < t < 12500 (fifth row, interval 
V). Note that y axis has not the same size in all panels. Right column: intersection 
points of orbit A with the PSS x — 0, p x > for the same time intervals (black points). 
In order to get a clear picture of the structural evolution of the phase space, in each 
panel the PSS corresponding to the central time of each time interval (t = 1250, 3750, 
6250, 8750 and 11250, from top to bottom) is plotted in gray. The orbit is initially 
regular and drifts from one island of stability to another, until finally its dynamical 
nature can be characterized as chaotic. 
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Figure 2. Time evolution of the logarithm of (a) the finite time MLE <7i, and (b) the 
reinitialized GALI 2 of orbit A. The five different time intervals I, II, III, IV and V, 
that correspond to the rows of figure [TJ are located between the vertical dashed gray 
lines. 



island of stability of figure DJb) and oscillates symmetrically around the bar's major- 
shown in figure Ufa) for < t < 2500 (interval I). Then in figure [H(d) we observe 
a first drift from the original island of stability to a nearby one for 2500 < t < 5000 
(interval II). The morphology of the orbit (figure (He)) changes at the same time to a 
different shape. 

In figures [D(e),(f) this transition has fully taken place and now the motion occurs 
on an island different from the one it started on, but its regular nature is still preserved 
for 5000 < t < 7500 (interval III). In the next time interval, however, 7500 < t < 10000 
(interval IV), we see a radical change in the orbit's morphology, which indicates the 
transition from regularity to chaoticity (figures CE^g),(h)). From figures |TJi),(j) we deduce 
that for 10000 <t< 12500 (interval V), orbit A moves to regions away from the bar in 
the configuration space, enters the big chaotic sea on the PSS, and from then on shows 
no regular behavior, as it remains in this chaotic region for the rest of the integration 
time. 

In figure [2(a) we depict the time evolution of the finite time MLE o\ (t) ( TTTT) of 
orbit A. During the first parts of the motion (intervals I and II), we clearly see a decay 
of 0i to zero indicating the regular nature of the orbit, in accordance with the results of 
figure [TJ Later on, at the end of interval III and mainly during interval IV, we witness 
a transient behavior where <j\ stops decaying and chaos arises. Then, in interval V <j\ 
remains positive and shows a tendency to slightly increase, which clearly suggests that 
orbit A becomes more chaotic as the bar's mass increases. 

Let us now examine the behavior of GALI2 for the same orbit. From figure E(b) 
we see that GALI 2 oscillates around a positive value during the time intervals I and II 
for which the orbit is regular and decays exponentially to zero, becoming < 10~ 8 , as 
soon as the orbit enters interval IV. In order to monitor the dynamical changes of the 
orbit we reinitialize the computation of GALI 2 as soon as GALI 2 < 10~ 8 , and plot in 
figure |3] the time needed for GALI 2 to decrease from GALI 2 = 1 to values smaller 
than 10~ 8 along the orbit's evolution. This figure demonstrates that orbit A is initially 
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Figure 3. The time intervals f<j needed for the reinitialized GALI2 to decrease from 
GALI2 = 1 to GALI2 < 10 W as a function of the integration time t of orbit A. The 
fluctuations of the td values (see insert) reflect the inhomogeneity of the dynamics in 
the corresponding time intervals. 



regular and its GALI2 becomes < 1CW for the first time after t ~ 9500 within interval 
IV. From that point on the orbit remains chaotic as its reinitialized GALI 2 repeatedly 
falls to zero very fast, resulting in small td values (td < 1000). This phase corresponds 
to the times that the orbit wanders in the big chaotic sea of the PSS (see figures QJh) 
and (j)). 

The important observation here is that the initial transition of orbit A from 
regularity to chaoticity and the subsequent variations of its dynamics are not distinctly 
captured by the evolution of the finite time MLE. Note from ( TlDl) and ( TTTT) that the 
MLE represents a time- averaged quantity over the whole evolution of the orbit, and 
consequently cannot reveal detailed changes of the orbit's motion as it exits an island 
and wanders within a large chaotic sea. On the other hand, GALI 2 does reveal such 
changes as it depends only on the current state of the dynamics and not on the previous 
history of the orbit. In fact, these advantages of the GALIs in capturing even brief 
dynamical transitions become more evident in the following sections, where we study 
orbits in 3 dof TD models . 

4-2. Orbits of the 3 dof model 

After describing how successful the finite time MLE and the GALI are in detecting 
changes in the chaotic vs. regular nature of orbits in the 2 dof restriction of Hamiltonian 
dl}, let us study some representative cases of the full 3 dof problem. 

4-2.1. A case where the bar gets weaker in time We now suppose that the bar's mass 
decreases linearly in time from an initial value Mb (to = 0) = 0.1 to Ms(t final — 
20000) = 0.0, following the law (jSJ) with a = —5 x 10~ 6 and study an orbit with quite 
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Figure 4. Projections of orbit B on the (a;, y) (left column) and the (x, z) plane 
(right column) for different time intervals: < t < 2500 (upper row, interval I), 
10000 < t < 12500 (middle row, interval II) and 17500 < t < 20000 (lower row, 
interval III). 



interesting behavior, which we call orbit B. Its initial condition is (x,y, z,p x ,p y ,p z ) = 
(0.3124,0.0,0.25,0.0,0.0), and its initial energy H B {t ) = E B (t ) = -0.429. In this 
case, the bar's strength starts at Qb{to) = 0.425 and reaches the value Qb{tfinai) = 
when there is no mass left at the bar component of the model. 

In figure H] we show the projections on the (x,y) (left column) and the (x,z) plane 
(right column) of orbit B for three different time intervals. Even by mere inspection one 
can observe the complexity of the evolution of orbit B. Although it is not safe to make 
accurate predictions for the nature of an orbit based on its form in the configuration 
space, we can say that orbit B looks regular in intervals I and III (upper and bottom 
rows of figure H] respectively), while it appears more complicated in interval II (middle 
row of figure H]). These observations suggest that the orbit is initially regular and after 
a chaotic phase becomes regular again. 

Exactly because orbit projections (especially for systems of more than 2 dof) are 
so difficult to interpret, we compute the o~\ and the GALI3 of orbit B, in order to 
analyze the stages through which the orbit passes. In figure 0(a) we see that <j\ 
decays for t < 5000, implying that the orbit is regular, then increases to higher values, 
indicating the possibility of chaotic motion, and finally for t > 14000 decays again to 
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Figure 5. Time evolution of the logarithm of (a) the finite time MLE ai, and (b) the 
reinitialized GALI 3 of orbit B. The three different time intervals I, II, and III, that 
correspond to the rows of figure 21 are located between the vertical dashed gray lines. 

zero, suggesting a return to regularity. 

It is remarkable how clearly GALI3 captures the different dynamical phases of the 
orbit's evolution (recall that that whenever GALI3 < 10~ 8 , we set GALI3 = 1 and repeat 
its computation using three orthonormal deviation vectors). GALI 3 initially oscillates 
around a high non-zero value, asserting that the motion is regular until about t ~ 5000 
when it starts to decrease rapidly to zero indicating that the orbit is chaotic. This phase 
lasts until t w 15000, when GALI3 jumps and remains practically constant until the end 
of the integration period, indicating that the orbit has become regular again. 

By comparing the panels of figure [5] we see that <j\ does not convincingly identify 
the transitions from regularity to chaoticity and vice versa, for reasons that were already 
discussed. Furthermore, we see that, even in the last interval of the orbit's evolution 
(t > 14000), where <j\ begins to fall towards zero, it decreases so slowly that the nature 
of the orbit is far from clear. On the other hand, GALI 3 can successfully detect the 
local properties of the dynamics and provide us with a clear knowledge of the chaotic 
vs. regular nature of an orbit, even for small time windows, where the orbit's nature 
often changes rapidly'. 

The distinction between the regular and chaotic intervals of orbit B is well depicted 
in figure |6l where the time t& needed for the reinitialized GALI 3 to become < 10 -8 is 
plotted as a function of the integration time. Small td values for 7500 < t < 14000 
correspond to chaotic epochs, where GALI 3 goes to zero exponentially fast, while larger 
td values correspond to intervals where GALI 3 takes longer to decay to zero. Observe 
also in figure [6] the remarkable fact (shown by an upwardly pointing arrow) that, after 
t > 15000, GALI3 no longer falls to zero (see figure [5]^b)) until the end of the integration 
time! This is certainly not expected in TI systems. It does occur, however, for orbit B 
of the 3 dof TD model, as well as orbit C of a similar system (see below). 

4-2.2. A case where the bar gets stronger in time Let us now study the case of a 
linear increase of the bar's mass Mb from the initial value Ms(to = 0) = 0.1 to 



Chaotic and regular motion in a time- dependent barred galaxy model 



14 



8000 
7000 
6000 
5000 
^ 4000 
3000 
2000 
1000 


2500 5000 7500 10000 12500 15000 

t 

Figure 6. The time td that the reinitialized GALI3 needs to decrease from GALI3 = 1 
to GALI3 < 10~ 8 , as a function of the integration time t of orbit B. Note that for 
t > 15000 GALI3 no longer falls to zero, indicating that the motion has entered a 
regular domain. 

Ms(t final = 20000) = 0.2, as we did in the case of the 2 dof model. We take again 
a = 5 x 10~ 6 in OH]). As an example, we consider the orbit C with initial condition 
(x,y,z, Px , Py ,p z ) = (0.225,0.0,0.25,0.0,0.0) and initial energy H c (t Q ) = E c (t ) = 
—0.441, which undergoes an interesting sequence of dynamical transitions. In this case, 
Qb starts from Qb(t ) = 0.425 and reaches the value Qb{t final) = 0.6732 at the end of 
the orbit's evolution. 

In figure [7] we plot the projections of orbit C on the (x, y) and the (x, z) planes 
(left and right column respectively) for three different time intervals. The orbit looks 
more or less regular in intervals I and III, although its shape is quite different in the 
two intervals. In interval II it looks a bit more complicated and seems to represent a 
transition between the two different configurations of intervals I and III. However, as 
has already been mentioned, the mere inspection of the orbit is not enough to accurately 
inform us about its chaotic or regular nature. 

In figure |8]^a) we plot the time evolution of 0\ for orbit C. From this figure we see 
that <j\ initially decays to zero, suggesting the regular character of the orbit. Then, 
at t ~ 5500 we observe a small increase of a\, which indicates a dynamical change in 
the orbit's behavior. This is soon followed by a monotonic decrease of o~i, which might 
indicate that orbit C becomes regular again. All this, however, is highly speculative. 
By contrast, the time evolution of the reinitialized GALI 3 (see figure E^b)) shows a 
lot more clearly the transition epoch between the two different regular states. Initially 
GALI3 remains constant and different from zero, providing clear evidence that the orbit 
is regular. Then, for 5000 < t < 8.500 it decreases to < 10~ 8 over a relatively long 
time interval, indicating a fundamental change in the character of the orbit. Finally, 
after reinitializing the index's computation at t « 8.500, GALI 3 begins to converges to 
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Figure 7. Projections of orbit C on the (x,y) (left column) and the (x, z) plane 
(right column) for different time intervals: < t < 2500 (upper row, interval I), 
5000 < t < 7500 (middle row, interval II) and 17500 < t < 20000 (lower row, interval 
III). 
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a positive constant, demonstrating the remarkable fact that the orbit has again become 
regular! Thus, here also, as in the case of orbit B (see figure G^b)), an interlude of 
chaotic behavior is followed by a transition to regularity, which lasts until the end of 
our integration time. 

4-3. Global dynamics of the 3 dof model 

After establishing the efficiency of the GALI method in identifying time intervals where 
an orbit of a TD model is regular or chaotic, let us use it to study in a more global 
way the dynamics of Hamiltonian (CQ). In particular, we will investigate the case 
considered in section 14.2.2] where the mass of the bar component increases linearly 
from M B (t = 0) = 0.1 to M B (t final = 20000) = 0.2, corresponding to a = 5 x 10~ 6 in 
©. 

In [37] the TI version of model ([1]) was considered for fixed values Mg = 0.1 
and Mb = 0.2 respectively. In that work, these two cases were referred as models 'S' 
and 'M' respectively. Ensembles of 50000 different initial conditions were integrated 
up to t = 10000 time units, and the GALI method was used to accurately determine 
the percentages of chaotic orbits. The analysis performed in [37J showed that chaotic 
behavior is more dominant for the 'M' model (i.e. the one with the more massive bar 
component). Our TD model (CQ) coincides with model 'S' of |37j at t — and becomes 
model 'M' at t = 20000. Thus, it is of interest to check for this model if the same sets of 
initial conditions considered in [37] show a tendency to increase their chaoticity at time 
grows from t — to t — 20000, in agreement with the general trend found in [37]. For 
this reason, we evolve the same three classes of initial condition distributions considered 
in [37J: 

• distribution /: 5000 orbits equally spaced in the space (x,z,p y ) with x G [0.0,7.0], 
z G [0.0, 1-5], p„ G [0.0,0.45] and (y,p x ,p z ) = (0,0,0), 

• distribution II: 5000 orbits equally spaced in the space (x,p y ,p z ) with x G [0.0, 7.0], 
p y G [0.0, 0.35], p z G [0.0,0.35] and (y,z, Px ) = (0,0,0), 

• distribution III: 5000 orbits whose spatial coordinates are chosen randomly over 
the mass density distribution of model 'S' (according to the so-called rejection 
method) within the rectangular box —a < x < a, —b < y < b, — c < z < c, 
with (p y ,Pz) = (0,0), and p x > obtained from (TjQ) with H taking for each initial 
condition a fixed random value in the interval [—0.22, 0] (see [37] for more details). 

The only difference between the distributions of our study and those of [37] is that we 
use only 5000 initial conditions, instead of 50000 used in [37], in order to facilitate our 
computations. 

Since our model is TD the percentage of chaotic orbits of the three orbital 
distributions is expected to change in time. In order to monitor these changes we adopt 
the following strategy: We divide the total integration time of 20000 time units in eight 
successive time windows of length At = 2500 time units. At the beginning of each 
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Figure 9. Time evolution of the percentages of chaotic orbits for the initial 
distributions I (dotted line), II (dashed line), III (solid line) of model ([T]), when 
M B is linearly increased from Mb {to = 0) = 0.1 to M B (tf lna i = 20000) = 0.2 (see 
text for more details). The percentages are obtained for successive time windows of 
length At = 2500 time units, and are attributed to the starting times of these intervals. 
For comparison, the percentages of the TI models 'S' (with Mb — 0.1) and 'M' (with 
Mb = 0.2) obtained in [37], are also plotted at t = and t = 20000 respectively, 
by filled squares (distribution /), filled circles (distribution II), and filled triangles 
(distribution III). 

time window, we reinitialize GALI3 to unity and follow a new set of three orthonormal 
deviation vectors for each orbit. Then for each time window we calculate the "current 
percentage of chaotic orbits" as the fraction of orbits whose GALI3 becomes < 10~ 8 
in the duration of the time window. The results of this procedure are shown in figure 
IHl where the percentages obtained for each time window are attributed to its starting 
time. The percentages found in [37] for models 'S' and 'M' are also plotted at t — and 
t = 20000 respectively, by filled squares (distribution J), filled circles (distribution i7), 
and filled triangles (distribution III). 

In accordance with the results of [37] we see an increase of the fraction of chaotic 
orbits for all three ensembles of orbits. Distribution III contains more chaotic orbits 
among the three ensembles, while the percentages of chaotic orbits is smallest for 
distribution 77, in agreement with what was observed in [37]. Obviously, one should 
not expect to find the same numerical values with the percentages reported in [37]. For 
example, our TD model is identical to model 'S' of [37] only momentarily for t — 0, 
and thus the percentages of the TD model reported for t = do not correspond exactly 
to model 'S'. Apart from the fact that the percentages reported in [37] were obtained 
for a TI model, another difference is that in [37] each orbit was integrated in a fixed 
Hamiltonian system for 10000 time units, while in our study each orbit is integrated for 
only 2500 time units in a TD potential. Nevertheless, the dynamical trends obtained 
in our study are in good agreement with the ones presented in [37], clearly showing the 



Chaotic and regular motion in a time- dependent barred galaxy model 



18 



efficiency of the GALI method in analyzing TD models. 
5. Conclusions 

Autonomous Hamiltonian systems are conveniently studied for fixed values of the total 
energy, where the location and extent of their regular and chaotic regions are time 
independent and can be accurately identified by a variety of methods especially in 
the low degree of freedom case. Even in such TI systems the dynamics can exhibit 
remarkable complexity, as there exist regimes of "strong" and "weak" chaos, as well 
as varying degrees of regularity, as the motion can occur on invariant tori of different 
dimensions and exhibit surprising localization properties in configuration and frequency 
space |12] . 

Naturally, therefore, Hamiltonian systems which are explicitly TD are expected to 
be a lot more complicated, since, in the absence of TI integrals, all the above attributes 
evolve in time. For example, in the TI case, orbits do not change their nature: If they 
are initially regular they will always remain so, while if chaotic they can get trapped 
for long times on the boundary of regular regimes (exhibiting a "weak" form of chaos), 
but will never entirely relinquish their chaoticity. This is not so in the TD case, where 
individual trajectories may indeed display sudden transitions from regular to chaotic 
behavior and vice versa during their time evolution. 

In the present paper, we have sought to shed some new light on these fascinating 
phenomena by studying the dynamics of a mean field model of a barred galaxy, whose 
mass parameters are allowed to vary linearly in time. Our primary goal was to show 
that transitions from order to chaos and vice versa do occur in this model and can be 
monitored much more accurately by local methods such as the GALI spectrum, rather 
than the more traditional approach of LEs. In addition, we wanted to investigate some 
astronomical properties of this TD system as it does incorporate some of the features 
appearing in N-body simulations and TI analytic potentials. 

In this regard, we have chosen for simplicity to vary only two parameters in time, 
Mb and Mp, keeping the size of the bar and the pattern speed fixed. Since the total 
mass is constant, whatever mass the bar loses is gained by the disc component and vice- 
versa. Moreover, to investigate more thoroughly the observed dynamical transitions, we 
have extended the maximal integration time of the orbits to T = 20000 Myr (20 billion 
years), which corresponds to a total interval of nearly 2 Hubble times. However, our 
results demonstrate that rich behavior and fundamental changes can also be observed 
within a single Hubble time of 10000 Myr. 

Most importantly, the GALI method turns to be again very efficient and accurate 
in the detection of chaotic motion in a TD system, as in the case of TI models. 
Furthermore, our work reveals that the method is especially suited for detecting intervals 
where an orbit changes its state fundamentally. By following the times that the GALIs 
require to fall to zero, one can describe in detail the orbit's successive passages from 
order to chaos and vice versa. 
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Finally, we focus on a more global astronomical study, where one would like to 
estimate qualitatively and quantitatively the relative fraction of regular and chaotic 
motion in such galaxy models. To this end, we choose different sets of initial condition, 
launch them in phase space and classify regular and chaotic orbits, depending on whether 
the GALI fluctuates around a non-zero value or falls exponentially to zero. We have 
thus been able to verify that the conclusion of an earlier publication on the TI model 
[37] that the percentage of chaos grows as the mass of the bar increases holds true in 
the TD case as well. It would be highly interesting to investigate these questions in 
more realistic models, where besides the mass parameters the rotation frequency of the 
galaxy is also allowed to vary accordingly. 

In closing, it is important to point out the advantages of the GALI method over 
the computation of the finite time MLE, as the GALIs do succeed in clearly capturing 
the transitions between different regular states and identifying the intermediate chaotic 
phases. By contrast, the manifestation of these different dynamical behaviors is 
much less pronounced in the time evolution of the MLE described by <j\ in this 
paper. Evidently, the practice of averaging Lyapunov exponents over an orbit's history 
smoothens out their fluctuations over short-lived events and gives them meaning only 
in the sense of the long time limit. 
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